Alpha diversity
Paths and libraries setting
# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")
# Load supplementary packages
packages <- c("multcomp", "openxlsx", "kableExtra", "dplyr")
invisible(lapply(packages, require, character.only = TRUE))Load phyloseq object after decontam
ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")Alpha diversity plots
Whole
# select whole samples
ps.whole <- subset_samples(ps.filter, Organ=="Whole")
ps.whole@sam_data %>% data.frame() -> test
# estimate alpha diversity
measures <- c("Observed", "Chao1", "Shannon")
p1 <- plot_richness(ps.whole,
x="Sample",
color="Strain",
measures=measures,
nrow = 1)
# extract data from richness plot to custom it
df <- p1$data
# change levels and order of species and location for plot
levels(df$Species) <- c("AA", "CP", "CQ")
df$Species <- factor(df$Species, levels = c("AA", "CQ", "CP"))
new_location <- c("Field - Guadeloupe",
"Laboratory - Slab TC (Wolbachia -)",
"Laboratory - Lavar",
"Field - Bosc",
"Field - Camping Europe")
df$Strain <- factor(df$Strain, levels = new_location)
levels(df$Strain)## [1] "Field - Guadeloupe" "Laboratory - Slab TC (Wolbachia -)"
## [3] "Laboratory - Lavar" "Field - Bosc"
## [5] "Field - Camping Europe"
# convert df into data.table and set min and max for y-axis
dat <- data.table::data.table(df)
dat[, y_min := value*0.1, by = variable]
dat[, y_max := value*1.1 , by = variable]
dat[dat$variable=="Observed", "y_min"] <- 10
dat[dat$variable=="Chao1", "y_min"] <- 10
dat[dat$variable=="Observed", "y_max"] <- 55
dat[dat$variable=="Chao1", "y_max"] <- 55
#dat[dat$variable=="Shannon", "y_max"] <- 3
# plot
p <- ggplot(dat,aes(Species,value, color=Strain)) +
facet_wrap(~ variable, drop=T, scale="free_y")+
geom_boxplot(outlier.colour = NA, alpha=0.8, position = position_dodge2(width=1, preserve="single"))+
geom_point(size=1, aes(shape=dat$Strain), position=position_dodge(width=0.7, preserve='total'))+
scale_shape_manual("Strain", values = c(16,17,17,16,16),
labels = c("Field - Guadeloupe",
expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), " -)")),
"Laboratory - Lavar",
"Field - Bosc",
"Field - Camping Europe"))+
scale_color_manual("Strain",
values=c("#00BF7D", "#5B6BF7", "#00B0F6", "#A3A500", "#F8766D"),
labels = c("Field - Guadeloupe",
expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), " -)")),
"Laboratory - Lavar",
"Field - Bosc",
"Field - Camping Europe"))+
labs(x="Species", y = "Alpha Diversity Measure")+
theme(legend.text.align = 0)
p # final plot with more space on y-axis
p_bis <- p +
geom_blank(aes(y = y_min)) +
geom_blank(aes(y = y_max))
p_bisWhole with rarefaction
# Rarefaction with sample.size=100
ps.whole2 <- ps.whole %>% rarefy_even_depth(rngseed=1, sample.size = 100)## `set.seed(1)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(1); .Random.seed` for the full vector
## ...
## 1 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## NP34
## ...
## 4OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
ps.whole2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 63 taxa and 66 samples ]
## sample_data() Sample Data: [ 66 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 63 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 63 reference sequences ]
ps.whole2@sam_data %>% data.frame() -> test2
compute_read_counts(ps.whole2)## [1] 6600
p2 <- alpha_div(ps.whole2, measures)
p2## Rarefaction with sample.size=1000
ps.whole3 <- ps.whole %>% rarefy_even_depth(rngseed=2, sample.size = 1000)## `set.seed(2)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(2); .Random.seed` for the full vector
## ...
## 4 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## NP29NP30NP34NP36
## ...
ps.whole3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 67 taxa and 63 samples ]
## sample_data() Sample Data: [ 63 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 67 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 67 reference sequences ]
ps.whole3@sam_data %>% data.frame() -> test3
compute_read_counts(ps.whole3)## [1] 63000
p3 <- alpha_div(ps.whole3, measures)
p3Ovaries
# select ovary samples
ps.ovary <- subset_samples(ps.filter, Organ=="Ovary")
sample_data(ps.ovary)$Read_depth <- sample_sums(ps.ovary)
# change factors and order in metadata
metadata.ps <- data.frame(sample_data(ps.ovary))
levels(metadata.ps$Strain)## [1] "Field - Bosc" "Field - Camping Europe" "Field - Guadeloupe"
## [4] "Laboratory - Lavar"
levels(metadata.ps$Strain) <- c(levels(metadata.ps$Strain), "Laboratory - Lavar")
metadata.ps$Strain[metadata.ps$Strain=="Lavar (labo)"] <- "Laboratory - Lavar"
levels(metadata.ps$Strain)## [1] "Field - Bosc" "Field - Camping Europe" "Field - Guadeloupe"
## [4] "Laboratory - Lavar"
metadata.ps$Strain <- droplevels(metadata.ps$Strain)
levels(metadata.ps$Strain)## [1] "Field - Bosc" "Field - Camping Europe" "Field - Guadeloupe"
## [4] "Laboratory - Lavar"
# update metadata in phyloseq object
sample_data(ps.ovary) <- metadata.ps
# estimate alpha diversity
p1 <- plot_richness(ps.ovary,
x="Sample",
color="Strain",
measures=measures,
nrow = 1)
# extract data from richness plot to custom it
df <- p1$data
# changer levels and order of species and location for plot
levels(df$Species) <- c("AA", "CP", "CQ")
df$Species <- factor(df$Species, levels = c("AA", "CQ", "CP"))
new_location <- c("Field - Guadeloupe",
"Laboratory - Slab TC (Wolbachia -)",
"Field - Bosc",
"Field - Camping Europe",
"Laboratory - Lavar")
df$Strain <- factor(df$Strain, levels = new_location)
# final plot
p_ovary <- ggplot(df,aes(Species,value,colour=Strain)) +
facet_wrap(~ variable, drop=T,scale="free")+
scale_color_manual("Strain",
values=c("#00BF7D", "#A3A500", "#F8766D", "#00B0F6"),
labels = c("Field - Guadeloupe",
"Field - Bosc",
"Field - Camping Europe",
"Laboratory - Lavar"))+
geom_boxplot(outlier.colour = NA, alpha=0.8, position = position_dodge2(width=1, preserve="single")) +
geom_point(size=2, position=position_dodge(width=0.7, preserve='total'))+
labs(x="Species", y = "Alpha Diversity Measure")+
theme(legend.text.align = 0)
p_ovaryMean of estimated alpha diversity
df_score <- data.frame(p$data)
x <- c("Culex pipiens - Field", "Culex pipiens - Bosc", "Culex pipiens - Camping Europe", "Culex pipiens - Lavar (lab)",
"Culex quinquefasciatus - Guadeloupe", "Culex quinquefasciatus - Slab TC",
"Aedes aegypti - Guadeloupe")
y <- c("Observed", "Chao1", "Shannon")
df <- data.frame(matrix(ncol=3, nrow=7))
rownames(df) <- x
colnames(df) <- y
# Aedes aegypti (from Guadeloupe)
df_score[df_score$Species=="AA" & df_score$variable=="Observed", "value"] %>% mean() -> df[7,1]
df_score[df_score$Species=="AA" & df_score$variable=="Chao1", "value"] %>% mean() -> df[7,2]
df_score[df_score$Species=="AA" & df_score$variable=="Shannon", "value"] %>% mean() -> df[7,3]
# Culex pipiens from field
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Field=="Field", "value"] %>% mean() -> df[1,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Field=="Field", "value"] %>% mean() -> df[1,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Field=="Field", "value"] %>% mean() -> df[1,3]
# Culex pipiens from labo
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,3]
# Culex pipiens from Bosc
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,3]
# Culex pipiens from Camping Europe
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,3]
# Culex quinquefasciatus from field (Guadeloupe)
df_score[df_score$Species=="CQ" & df_score$variable=="Observed" & df_score$Field=="Field", "value"] %>% mean() -> df[5,1]
df_score[df_score$Species=="CQ" & df_score$variable=="Chao1" & df_score$Field=="Field", "value"] %>% mean() -> df[5,2]
df_score[df_score$Species=="CQ" & df_score$variable=="Shannon" & df_score$Field=="Field", "value"] %>% mean() -> df[5,3]
# Culex quinquefasciatus from labo (Wolbachia-)
df_score[df_score$Species=="CQ" & df_score$variable=="Observed" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,1]
df_score[df_score$Species=="CQ" & df_score$variable=="Chao1" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,2]
df_score[df_score$Species=="CQ" & df_score$variable=="Shannon" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,3]
df %>%
kbl(booktable=TRUE) %>%
kable_paper("hover", full_width = F)| Observed | Chao1 | Shannon | |
|---|---|---|---|
| Culex pipiens - Field | 20.84211 | 22.17982 | 0.1971403 |
| Culex pipiens - Bosc | 20.61538 | 20.94231 | 0.1790446 |
| Culex pipiens - Camping Europe | 21.33333 | 24.86111 | 0.2363475 |
| Culex pipiens - Lavar (lab) | 27.90909 | 30.03485 | 1.0724543 |
| Culex quinquefasciatus - Guadeloupe | 17.25000 | 24.66667 | 1.5301268 |
| Culex quinquefasciatus - Slab TC | 28.46154 | 31.08269 | 1.8549403 |
| Aedes aegypti - Guadeloupe | 32.22222 | 33.51852 | 1.0625983 |
Save plots
tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.tiff", units="in", width=8, height=5, res=300)
p
dev.off()## quartz_off_screen
## 2
png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.png", units="in", width=8, height=5, res=300)
p
dev.off()## quartz_off_screen
## 2
pdf ( "../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.pdf")
p_bis
dev.off()## quartz_off_screen
## 2
tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_100.tiff", units="in", width=8, height=5, res=300)
p2
dev.off()## quartz_off_screen
## 2
png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_100.png", units="in", width=8, height=5, res=300)
p2
dev.off()## quartz_off_screen
## 2
tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_1000.tiff", units="in", width=8, height=5, res=300)
p3
dev.off()## quartz_off_screen
## 2
png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_1000.png", units="in", width=8, height=5, res=300)
p3
dev.off()## quartz_off_screen
## 2
Statistics
Prepare subdatasets for stats
# Whole
df <- p$data
df_culex <- df[df$Species!="AA",]
df_pipiens <- df[df$Species=="CP",]
df_quinque <- df[df$Species=="CQ",]
# Whole rarefied (100)
df2 <- p2$data
df2_culex <- df2[df2$Species!="AA",]
df2_pipiens <- df2[df2$Species=="CP",]
df2_quinque <- df2[df2$Species=="CQ",]
# Whole rarefied (1000)
df3 <- p3$data
df3_culex <- df3[df3$Species!="AA",]
df3_pipiens <- df3[df3$Species=="CP",]
df3_quinque <- df3[df3$Species=="CQ",]
# Wolbachia- vs Wolbachia+
df_culex_wolbachia <- df_culex[df_culex$Species=="CQ" & df_culex$Strain!="Laboratory - Lavar",]
# Wolbachia- vs Wolbachia+ (rarefied 100)
df2_culex_wolbachia <- df2_culex[df_culex$Species=="CQ" & df2_culex$Strain!="Laboratory - Lavar",]
# Wolbachia- vs Wolbachia+ (rarefied 1000)
df3_culex_wolbachia <- df3_culex[df_culex$Species=="CQ" & df3_culex$Strain!="Laboratory - Lavar",]Influence of depth sequencing
Effect of Species
Observed
# names for final array
names <- c("Group", "Read_depth", "Group:Read_depth", "Residuals")
names2 <- c("Group", "Read_depth", "Residuals")
## Observed
df_observed <- df[df$variable=="Observed" & df$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
wb <- createWorkbook()
addWorksheet(wb, "Alpha diversity (depth seq)")
writeData(wb, "Alpha diversity (depth seq)", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 6, rowNames = TRUE)
## Observed (rarefied 100)
df_observed_100 <- df2[df$variable=="Observed" & df2$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 16, rowNames = TRUE)
## Observed (rarefied, 1000)
df_observed_1000 <- df3[df3$variable=="Observed" & df3$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 26, rowNames = TRUE)Chao1
# Chao1
df_chao <- df[df$variable=="Chao1" & df$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 6, rowNames = TRUE)
# Chao1 (rarefied, 100)
df_chao_100 <- df2[df2$variable=="Chao1" & df2$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 16, rowNames = TRUE)
## Chao (rarefied, 1000)
df_chao_1000 <- df3[df3$variable=="Chao1" & df3$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 1000)", startCol = 9, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 26, rowNames = TRUE)
#saveWorkbook(wb, "Supplementary_Table_2_review_sheet1_30_11_21.xlsx", overwrite = TRUE)Shannon
## Shannon
df_shannon <- df[df$variable=="Shannon" & df$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 6, rowNames = TRUE)
## Shannon (rarefied, 100)
df_shannon_100 <- df2[df2$variable=="Shannon" & df2$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 16, rowNames = TRUE)
## Shannon (rarefied, 1000)
df_shannon_1000 <- df3[df3$variable=="Shannon" & df3$Field!="Lab ",]
#### Construct linear model
model <- lm(value~ Species*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 26, rowNames = TRUE)
saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of tetracycline
Observed
## Observed
df_observed<- df_quinque[df_quinque$variable=="Observed",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 35)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 40, rowNames = TRUE)
## Observed (rarefied 100)
df_observed_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Observed",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 50, rowNames = TRUE)
## Observed (rarefied, 1000)
df_observed_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Observed",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 60, rowNames = TRUE)Chao1
## Chao1
df_chao <- df_quinque[df_quinque$variable=="Chao1",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 40, rowNames = TRUE)
## Chao1 (rarefied 100)
df_chao_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Chao1",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 50, rowNames = TRUE)
## Chao1 (rarefied, 1000)
df_chao_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Chao1",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
#ggpubr::ggqqplot(residuals(model))
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO INDEX (rarefied, 1000)", startCol = 9, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 60, rowNames = TRUE)Shannon
## Shannon
df_shannon <- df_quinque[df_quinque$variable=="Shannon",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
#ggpubr::ggqqplot(residuals(model))
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 40, rowNames = TRUE)
## Shannon (rarefied 100)
df_shannon_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Shannon",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 50, rowNames = TRUE)
## Shannon (rarefied, 1000)
df_shannon_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Shannon",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 60, rowNames = TRUE)
saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)Effect of laboratory
Observed
## Observed
df_observed <- df_pipiens[df_pipiens$variable=="Observed",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "Effect of lab and location in the Culex pipiens samples", startCol = 1, startRow = 68)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 72, rowNames = TRUE)
## Observed (rarefied 100)
df_observed_100 <- df2_pipiens[df2_pipiens$variable=="Observed",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 82, rowNames = TRUE)
## Observed (rarefied, 1000)
df_observed_1000 <- df3_pipiens[df3_pipiens$variable=="Observed", ]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 92, rowNames = TRUE)Chao1
## Chao1
df_chao <- df_pipiens[df_pipiens$variable=="Chao1",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 72, rowNames = TRUE)
## Chao1 (rarefied 100)
df_chao_100 <- df2_pipiens[df2_pipiens$variable=="Chao1",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 82, rowNames = TRUE)
## Chao1 (rarefied, 1000)
df_chao_1000 <- df3_pipiens[df3_pipiens$variable=="Chao1", ]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO INDEX (rarefied, 1000)", startCol = 9, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 92, rowNames = TRUE)Shannon
## Shannon
df_shannon <- df_pipiens[df_pipiens$variable=="Shannon",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 72, rowNames = TRUE)
## Shannon (rarefied 100)
df_shannon_100 <- df2_pipiens[df2_pipiens$variable=="Shannon",]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 82, rowNames = TRUE)
## Shannon (rarefied, 1000)
df_shannon_1000 <- df3_pipiens[df3_pipiens$variable=="Shannon", ]
#### Construct linear model
model <- lm(value~ Strain*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")
#### Create residuals plot
par(mfrow=c(2,2))
plot(model)#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 92, rowNames = TRUE)
saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)Analysis of alpha diversity
Effect of species
Observed
names <- c("Species", "Residuals")
## Observed
df_observed <- df[df$variable=="Observed" & df$Field!="Lab ",]
cat("OBSERVED\n\n")## OBSERVED
res <- make_anova_tukey(df_observed, df_observed$Species, df_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 974.89 487.44 10.543 0.0003622 ***
## Residuals 29 1340.83 46.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## CQ - AA == 0 -14.972 4.086 -3.664 0.00265 **
## CP - AA == 0 -11.380 2.751 -4.136 < 0.001 ***
## CP - CQ == 0 3.592 3.741 0.960 0.60166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
addWorksheet(wb, "Alpha diversity")
writeData(wb, "Alpha diversity", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 12, rowNames = TRUE)Chao1
## Chao1
df_chao <- df[df$variable=="Chao1" & df$Field!="Lab ",]
cat("CHAO1\n\n")## CHAO1
res <- make_anova_tukey(df_chao, df_chao$Species, df_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 789.86 394.93 8.5122 0.001234 **
## Residuals 29 1345.47 46.40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## CQ - AA == 0 -8.852 4.093 -2.163 0.0923 .
## CP - AA == 0 -11.339 2.756 -4.114 <0.001 ***
## CP - CQ == 0 -2.487 3.747 -0.664 0.7824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 12, rowNames = TRUE)Shannon
## Shannon
df_shannon <- df[df$variable=="Shannon" & df$Field!="Lab ",]
cat("SHANNON\n\n")## SHANNON
res <- make_anova_tukey(df_shannon, df_shannon$Species, df_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 8.4685 4.2342 46.442 9.09e-10 ***
## Residuals 29 2.6440 0.0912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## CQ - AA == 0 0.4675 0.1814 2.577 0.0384 *
## CP - AA == 0 -0.8655 0.1222 -7.083 <0.001 ***
## CP - CQ == 0 -1.3330 0.1661 -8.025 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 12, rowNames = TRUE)Effect of tetracycline
names <- c("Strain", "Residuals")
## Observed
culex_observed <- df_quinque[df_quinque$variable=="Observed",]
cat("OBSERVED\n\n")## OBSERVED
res <- make_anova_tukey(culex_observed, culex_observed$Strain, culex_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 1 384.49 384.49 21.521 0.0003211 ***
## Residuals 15 267.98 17.87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 11.212
## Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 2.417 4.639
## Pr(>|t|)
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.000321 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 20)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 31, rowNames = TRUE)Chao1
## Chao1
culex_chao <- df_quinque[df_quinque$variable=="Chao1",]
cat("CHAO1\n\n")## CHAO1
res <- make_anova_tukey(culex_chao, culex_chao$Strain, culex_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 1 125.92 125.918 5.3198 0.03577 *
## Residuals 15 355.04 23.669
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 6.416
## Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 2.782 2.306
## Pr(>|t|)
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.0358 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 31, rowNames = TRUE)Shannon
## Shannon
culex_shannon <- df_quinque[df_quinque$variable=="Shannon",]
cat("SHANNON\n\n")## SHANNON
res <- make_anova_tukey(culex_shannon, culex_shannon$Strain, culex_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 1 0.32272 0.32272 1.7577 0.2047
## Residuals 15 2.75404 0.18360
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.3248
## Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.2450 1.326
## Pr(>|t|)
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.205
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 31, rowNames = TRUE)Effect of lab in Culex pipiens samples
Observed
## Observed
pipiens_observed <- df_pipiens[df_pipiens$variable=="Observed",]
cat("OBSERVED\n\n")## OBSERVED
res <- make_anova_tukey(pipiens_observed, pipiens_observed$Strain, pipiens_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 511.28 255.642 13.195 4.45e-05 ***
## Residuals 38 736.23 19.374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0 -7.2937 1.5398 -4.737
## Field - Camping Europe - Laboratory - Lavar == 0 -6.5758 2.0272 -3.244
## Field - Camping Europe - Field - Bosc == 0 0.7179 2.1724 0.330
## Pr(>|t|)
## Field - Bosc - Laboratory - Lavar == 0 < 0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0 0.00682 **
## Field - Camping Europe - Field - Bosc == 0 0.94077
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "Effect of lab (and location) in the Culex pipiens samples", startCol = 1, startRow = 37)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 47, rowNames = TRUE)Chao1
## Chao1
pipiens_chao <- df_pipiens[df_pipiens$variable=="Chao1",]
cat("CHAO1\n\n")## CHAO1
res <- make_anova_tukey(pipiens_chao, pipiens_chao$Strain, pipiens_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 692.1 346.05 10.453 0.0002415 ***
## Residuals 38 1258.0 33.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0 -9.093 2.013 -4.517
## Field - Camping Europe - Laboratory - Lavar == 0 -5.174 2.650 -1.952
## Field - Camping Europe - Field - Bosc == 0 3.919 2.840 1.380
## Pr(>|t|)
## Field - Bosc - Laboratory - Lavar == 0 <0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0 0.135
## Field - Camping Europe - Field - Bosc == 0 0.357
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 47, rowNames = TRUE)Shannon
## Shannon
pipiens_shannon <- df_pipiens[df_pipiens$variable=="Shannon",]
cat("SHANNON\n\n")## SHANNON
res <- make_anova_tukey(pipiens_shannon, pipiens_shannon$Strain, pipiens_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res## $anova
## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## var 2 7.8247 3.9124 18.859 2.047e-06 ***
## Residuals 38 7.8832 0.2075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $tukey
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = value ~ var, data = df)
##
## Linear Hypotheses:
## Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0 -0.8934 0.1593 -5.607
## Field - Camping Europe - Laboratory - Lavar == 0 -0.8361 0.2098 -3.986
## Field - Camping Europe - Field - Bosc == 0 0.0573 0.2248 0.255
## Pr(>|t|)
## Field - Bosc - Laboratory - Lavar == 0 <0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0 <0.001 ***
## Field - Camping Europe - Field - Bosc == 0 0.964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 47, rowNames = TRUE)
saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)Session info
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.3.1 openxlsx_4.2.3 multcomp_1.4-14 TH.data_1.0-10
## [5] MASS_7.3-53 survival_3.2-7 mvtnorm_1.1-1 forcats_0.5.0
## [9] stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4 readr_1.4.0
## [13] tidyr_1.1.2 tibble_3.0.4 tidyverse_1.3.0 ggplot2_3.3.2
## [17] phyloseq_1.30.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-149 fs_1.5.0 lubridate_1.7.9
## [4] webshot_0.5.2 httr_1.4.2 tools_3.6.3
## [7] backports_1.1.10 bslib_0.2.4 R6_2.4.1
## [10] vegan_2.5-6 DBI_1.1.0 BiocGenerics_0.32.0
## [13] mgcv_1.8-33 colorspace_2.0-0 permute_0.9-5
## [16] ade4_1.7-15 withr_2.3.0 tidyselect_1.1.0
## [19] compiler_3.6.3 cli_2.1.0 rvest_0.3.6
## [22] Biobase_2.46.0 xml2_1.3.2 sandwich_3.0-0
## [25] labeling_0.4.2 bookdown_0.22 sass_0.3.1
## [28] scales_1.1.1 digest_0.6.26 rmarkdown_2.7
## [31] XVector_0.26.0 pkgconfig_2.0.3 htmltools_0.5.1.1
## [34] highr_0.8 dbplyr_1.4.4 rlang_0.4.10
## [37] readxl_1.3.1 rstudioapi_0.11 farver_2.0.3
## [40] jquerylib_0.1.3 generics_0.0.2 zoo_1.8-8
## [43] jsonlite_1.7.1 zip_2.1.1 magrittr_1.5
## [46] biomformat_1.14.0 Matrix_1.2-18 fansi_0.4.1
## [49] Rcpp_1.0.5 munsell_0.5.0 S4Vectors_0.24.4
## [52] Rhdf5lib_1.8.0 ape_5.4-1 lifecycle_0.2.0
## [55] stringi_1.5.3 yaml_2.2.1 zlibbioc_1.32.0
## [58] rhdf5_2.30.1 plyr_1.8.6 grid_3.6.3
## [61] blob_1.2.1 parallel_3.6.3 crayon_1.3.4
## [64] lattice_0.20-41 Biostrings_2.54.0 haven_2.3.1
## [67] splines_3.6.3 multtest_2.42.0 hms_0.5.3
## [70] ps_1.4.0 knitr_1.30 pillar_1.4.6
## [73] igraph_1.2.6 reshape2_1.4.4 codetools_0.2-16
## [76] stats4_3.6.3 reprex_0.3.0 glue_1.4.2
## [79] evaluate_0.14 data.table_1.13.2 modelr_0.1.8
## [82] vctrs_0.3.4 rmdformats_1.0.2 foreach_1.5.1
## [85] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [88] xfun_0.22 broom_0.7.2 viridisLite_0.3.0
## [91] iterators_1.0.13 IRanges_2.20.2 cluster_2.1.0
## [94] ellipsis_0.3.1